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Recently I proposed [3| a new method for solving the pairing Hamiltonian with the pair-condensate 
wavefunction ansatz based on the Heisenberg equations of motion for the density matrix operators. 

In this work an improved version is given by deriving the relevant equations more carefully. I evaluate 
both versions in a large ensemble with random interactions, and the accuracy of the methods is 
given statistically in terms of root-mean-square derivations from the exact results. The widely used 
variational calculation is also done and the results and computing-time costs are compared. 
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I. INTRODUCTION II. FORMALISM 


Pairing correlation has long been recognized in nuclei 
0 and influences all the properties of the latter, such 
as mass, gap of excitation energy, and moment of in¬ 
ertia 0 . In general, any mean-field treatment of the 
nuclei needs to account it somehow to get reasonable re¬ 
sults. Among the methods the most popular one may 
be the Bardeen-Cooper-Schrieffer (BCS) theory i , or 
its advanced version the Hartree-Fock-Bogoliubov (HBF) 
theory Q , where pairing correlation is considered by in¬ 
troducing quasi-particles and writing the ground state as 
a vacuum of the latter. But there are disadvantages of 
breaking the exact particle number and need for a un¬ 
physical minimum pairing strength ill- 

A common improvement is to use the pair-condensate 
[Eq. dD), the BCS wavefunction projected onto good 
particle number] as the ground state wavefunction Q. 
Usually the criteria to determine the variational param¬ 
eters is minimizing the energy in the variation princ^le 
(variation after projection BCS) [Ml. In Ref. Q I 
proposed a new criteria based on the Heisenberg equa¬ 
tions of motion (EOM) for density matrix operators, as 
the lowest-order (mean-held) results of the generalized 
density matrix (GDM) formalism (l^ - [^ . The method 
was applied to the calcium isotopes with the well known 
FPD6 interaction [ 2 l|. 

In this work I give an improved version of the method. 
In Ref. [ij the relevant equations are derived assuming 
that neighboring even-even nuclei have the same Hartree- 
Fock (HF) single-particle energies and occupation num¬ 
bers; here I do it more carefully allowing the latters to be 
different. To see the validity of the approaches, both ver¬ 
sions are applied to a large ensemble with 1000 examples 
and random interactions; the results are good in almost 
all the examples. I also compare the results and nec¬ 
essary computing-time costs with those of the variation 
principle. In Sec. Hi derive the improved version of the 
GDM pairing theory. Then in Sec. IIIII the approaches 
are evaluated in the large random ensemble. Finally Sec. 
CYI summarizes the work and discusses further directions. 


The pairing theory in Ref. [T| was derived as the 
lowest-order (mean-held) results of the complete GDM 
formalism. In fact the derivation could be simpler and 
more clear if we focus on the mean helds and do not in¬ 
troduce collective (quadrupole) phonons. Also, in Ref. 
[l| I made the assumption that neighboring even-even 
nuclei have the same HF single-particle energies and oc¬ 
cupation numbers. Below I derive an improved version 
of the pairing theory abandoning the above assumption. 

As before, the ground state of the 2 iV-particle system 
is assumed to be an A-pair condensate, 

|(/.^) = ^(-P^)^|0), (1) 

\/Mv 

where xn is the normalization factor, and Pl is the pair 
creation operator 
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In Eq. ([ 2 ]) the summation runs over the entire single¬ 
particle space. |1) is the time-reversed level of the single¬ 
particle level | 1 ). The pair structure vi are parameters 
to be determined by the theory. 

With the antisymmetrized fermionic Hamiltonian 

H = ei 2 <i\a 2 + ^ X! 1^123401020304, ( 3 ) 

I calculate the exact Heisenberg equations of motion for 
the density matrix operators R12 = ojoi and Ari2 = 
0201 , 

[Ri 2 ,H] = [e,R]i 2 

--^ U 54320 ^o|a 3 ai + K X 1^1345020^0405, ( 4 ) 

^ 345 ^ 345 

[Ari2, H] = {eK)i 2 — {eK) 2 i + - Vi 234 Ar 43 

34 

- - ^ 145430^040302 + K X 1^25430^040301. ( 5 ) 

^ 345 ^ 345 
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Terms on the right-hand sides of Eqs. 0 and 0 are 
read as matrix multiplications; for example, [e,i?]i 2 = 
(ei?)i 2 — (i?e)i 2 = ei3^32 — Rl3£32- 

On the pair condensate o, the density matrices are 
“diagonal”: 

Pi2 = {4>N\alai\4>N) = Si2ni , ( 6 ) 

Ki 2 = 10201 |(/>Ar) = i 5 j 2 sf. (7) 

In practical shell-model calculations usually each single¬ 
particle level has distinct spin and parity, thus the mean 
fields are also “diagonal”: 

fl2 =£12 + ^ Vi432P^ = 6i2ei , (8) 

34 

^ • ( 9 ) 

34 

The pairing mean field should not be mixed with the 
Kronecker delta <5. 

Now I take matrix elements of the exact equa¬ 
tions of motion between the pairing ground states 
O: equation © between those of the same nu¬ 
clei (“(^Afl” and equation dS]) between 

those of neighboring even-even nuclei and 

“I^at)”). On the left-hand side of Eq. ([5]) 
we have {(I)n-i\[Ki 2 , H]\(j)]\r) = {(j)N_i\Ki 2 H\(j)N) - 
((l)N-l\HKi2\(j)N) ~ (En — En-i)k!^2^ where En and 
En -1 are the ground state energies, H\<j)N) ~ En\4>n) 
and H\(f>N-i) ~ En-i\(I)n-i)■ Similarly for Eq. (j4]) we 
have (()!)Ar|[i?i 2 ,iJ]|(/>Ar) ~ {En - En)p ^2 = 0 - the 
right-hand sides the “two-body density matrices” are ap¬ 
proximated in the following way: 

((/)Ar 104030201 |(/)7v) ~ (</'w|a 40 i|(/)Ar)((/'Ar| 0302 |(/)Ar) 

— ((/lAT 104021 (/) 7 v) ( 07 VI03O110jv) 
-|-(07v|oJa^|0Ar-i)(0Ar-i|o2ai|0Ar), 

= Pl4P^ - PmPi3 + «12: (10) 

(07V-1 10^030201 |07v) 

Ri (07V-l|a|oi|07V-l)(0Ar-l|o3O2|07v) 

— (0Ar-l|a4O2|0Ar-l)(07V-l|o3Oi|0Ar) 

+ (0Ar-l|a4O3|0Ar-l)(07V-l|o2Oi|0Ar) 

_ 7V-1 N N-1 N I N-1 N /I IN 

— P44 K23 P24 ^13 + 7*34 ^12- (11) 

Equations (nni and m would be exact if the ground 
states were taken as single-particle Slater determinants: 
the right-hand sides were just the fully contracted terms 
in Wick’s theorem. Here they are approximations be¬ 
cause the ground states are taken as pair condensates 
dH). Einally Eqs. dH and ([1]) become 

0=[fN,PN]- I^NSlf + SN^lf, (12) 

{En — En-i)kn = Jn-ikn + i^n/n-i 

+6n — SnPn-i — Pn-iSn- (13) 


Under the “diagonal” properties dfi]), ©, ©, and (O, Eq. 
dm is satisfied automatically, and Eq. m becomes 

Ejv — Ejv-i = 2e(^ ^ + Pi —iy-. (14) 

sf 

Equation (jl4p is the main equation of the improved the¬ 
ory. It implies that the right-hand side is independent of 
the single-particle label 1, by which the parameters vi in 
Eq. (2) are fixed. The main equation of the old theory, 
Eq. (19) in Ref. [l|, corresponds to replacing and 
by and rii . Equation m includes the well- 
known particle-particle random phase approximation Q 
as its special case of IV = 1 (n° = 0, e° = e). The nor¬ 
malization factor xn o, occupation numbers 0, 
pair-transfer amplitudes sf dll), mean fields ef dSl) ^’^d 
0 are functions of the pair structures vi 0; their 
functional forms, as the “kinematics” of the systeim have 
already been given in Eqs. (23) and (24) of Ref. [I| and 
are not repeated here. 

In the next section I take the pairing Hamiltonian: 

£12 = <5i2£1i ki234 = (15) 

in Eq. ©■ Consequently the mean fields 0 and 0 
become 

ef = £i-G'iinf, 9i =^^Gi 2 S 2 . (16) 
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I apply both versions of the theory to a large random 
ensemble; for convenience I call Eq. (fl^ the “GDM2” 
pairing theory, and Eq. (19) in Ref. [11 the “GDMl” 
pairing theory. 

III. RANDOM ENSEMBLE 

The “GDMl” theory was applied in Ref. [I| to calcium 
isotopes with the FPD6 interaction and the results are 
good. Here 1 would like to consider both the “GDMl” 
and “GDM2” theories in a large ensemble with random 
interactions; consequently we can speak statistically the 
accuracy of the theories in terms of the root-mean-square 
derivations from the exact results. The variational calcu¬ 
lation with the trial wavefunction o is also performed 
and the accuracy and computing-time cost are compared. 

The random ensemble has 1000 examples with differ¬ 
ent parameters determined in the following way. For 
each example, 1 first pick up the single-particle levels 
randomly from the pool 2j = 1,3,5,7,9,11,13. Each 
angular momentum j has a 40% probability of being 
selected; the selected ones (at least two) constitute the 
model space. Second, the single-particle energies ei m 
are determined randomly following the uniform distribu¬ 
tion from —10 MeV to 0. Third, 1 pick up the “pairing 
strength” Gmax as a random number following the uni¬ 
form distribution from 0 to 2 MeV, then the pairing ma¬ 
trix elements G 12 m are distributed uniformly from 0 
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to Gniax- Finally, the number of pairs N is determined 
following the uniform distribution from 1 to — 1, where 
2n = + 1) is the maximal particle number allowed 

by the model space. 

I perform four sets of calculations for the ensemble: 
two GDM calculations “GDMl” and “GDM2”, varia¬ 
tional calculation “VAR”, and the exact calculation. (In 
this work the exact calculation is done by diagonaliza- 
tion in spaces with fixed seniority [l^, [2^ . It can also be 
achieved by the Monte Carlo algorithm [2jj2g; or the 
Richardson’s method in some special cases [27H^ .l For 
the variational calculation of the pair-transfer amplitudes 
Si = (^Ar-i|ajai|0Ar) 0, I show two sets of results. In 
“VARl” the pair structure vi in \(j)N-i) and \(j)N) 
are the same, given by minimizing {(j)N\H\(f)]\[); while 
in “VAR2” vi in |(()Ar_i) and \(j)]sf) are different, given 
by minimizing {(j)N-i\H\(j)N-i) and {(j)N\H\4>N), respec¬ 
tively. 

In Figs. [T] and [5] I show the complete spectroscopic 
results for the ensemble. For example, in panel (a) of 
Fig. [T] there are 3154 points (crosses), corresponding to 
the 3154 single-particle levels in the 1000 examples of 
the ensemble. The horizontal coordinate of each point is 
the exact value of ni of the corresponding single-particle 
level, while the vertical coordinate is the “GDMl” value. 
Thus a perfect calculation would have all the points lying 
on the y = X straight line. Similarly, in panels (b) and (c) 
of Fig. [T]the vertical coordinates are the “GDM2” and 
“VAR” values of ni, respectively. Figure [His plotted in 
the same way for the pair-transfer amplitudes si. From 
Figs. |T]and|Hwe see that the variational calculation for 
ni and the “VAR2” version of si are generally better 
than the GDM ones. However, the less-careful “VARl” 
calculation of si is worse than the GDM ones. The root- 
mean-square (cr) derivations from the exact results are 

^GDMi ^ 0.0122, ^ 0.0188, = 0.0045,(17) 

on Fig. [T]and 

^ 0.0207, ^ 0.0191, 

^VARi ^ 0.0373, ^ 0.0091, (18) 

on Fig. IH We note that on Fig. |H there is a point 

much worse than others for all the calculations, so let 

us look at the particular example it belongs to. This 

example has 2N = 12 particles on two single-particle 

levels with angular momenta j = | and energies 

£5 = —8.988, £9 = —9.390 MeV. The pairing two-body 

matrix elements are Gb s = 0.4252, Ga a = 0.0456, and 

Gsa = Gaa = 0.0016 MeV. The failed point corre- 
2'2 2 ’ 2 

sponds to the j = f level. This example is particu¬ 
lar in that £5 — £a ~ Gs 5 ^ Ga 2 ^ Gs a, thus 
2 2 2’2 2’2 2’2 

there is little correlation between the two levels. I have 

looked at the exact wavefunction of the daughter nucleus 

and it is mainly Pj (Pg )"'^|0), which is not representable 
2 2 

by Eq. O- Without this example, the root-mean- 
square derivations for the pair-transfer amplitudes si are 


^GDMi ^ 0.0177, ^ 0.0163, = 0.0352, and 

^yAR2 ^ 0.0079. 

The necessary formula for the GDM and variational 
calculations were given in Ref. [l|. In general, in large 
model spaces the GDM calculation costs less time than 
the variational one, by a factor of the number of non¬ 
degenerate single-particle levels (in order of magnitude), 
because the former needs to calculate only {(j)N-i\Pi\(j)]\[) 
while the latter needs (^!)Ar|P/P 2 |(/>Ar) {Pi = aloj)- How¬ 
ever, as we can see from Figs. |T] and [Hand Eqs. (ITTl) and 
(ITSl) , the accuracy of the GDM method is close to that of 
the variational one. The reduction of time cost may be 
a big advantage when doing ab-initio mean-field calcula¬ 
tions for medium and heavy nuclei, especially if we were 
fitting parameters of the interaction (for example the ef¬ 
fort in developing density functionals with spectroscopic 
accuracy). 

Next I would like to see the accuracy of the meth¬ 
ods depending on different quantities. For this purpose 
I plot the root-mean-square derivations of ni and si as 
functions of the single-particle angular momentum j, par¬ 
ticle number 2iV, and pairing strength Gmax in Figs. [H 
m and |5l respectively. We have the following observa¬ 
tions: 1. From Figs. [3] and Sj we see that as a trend the 
GDM results improve with increasing j and iV, which 
should be expected because the GDM formalism is a col¬ 
lective theory. 2. There seems to be no obvious trend 
with the pairing strength shown in Fig. |5l 3. On 

Fig. m the GDM2 calculation for iV = 1, which is the 
particle-particle random phase approximation, is much 
better than the GDMl calculation. And for small N 
(until TV = 6), the GDM2 si seems to be slightly bet¬ 
ter than the GDMl si. 4. On Fig. [H the GDM2 si is 
slightly better than the GDMl si at the smallest angu¬ 
lar momenta j = \ and |. For these levels the difference 
Ui — , which is inversely proportional to j around 

the Fermi surface, is largest; and the GDM2 theory with 
seems to be slightly better than the GDMl one 
with Ui . 5. The GDMl ni seems to be consistently 
better than the GDM2 m and the reason is still unclear. 
6. In panel (c) of Figs. [H |H andjSjthe worst points are 
at 2j = 5, TV = 6, and 0.5 MeV < G^ax < 0.6 MeV, 
respectively; because these three sub-groups contain the 
“worst” example mentioned below Eq. m- 

At last in Fig. [51 1 show the results for the ground 
state energy by different calculations. In panel (a) there 
are 1000 points (crosses), corresponding to the 1000 ex¬ 
amples of the ensemble. The horizontal coordinate of 
each point is the pairing correlation energy of the cor¬ 
responding example, Apair = ~ PsheW, where 

Esheii is the exact ground state energy of the shell model 
calculation, and nf = 1 or 0 is the occupation num¬ 
ber of the naive Fermi distribution. The vertical coor¬ 
dinate shows the variational ground state energy mea¬ 
sured from the exact one, Evar = ~ Asheii, 

where is the pair condensate 0 with its pair 

structure vi 0 determined by the variation principle. 
Similarly for panels (b) and (c), but the vertical co- 
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ordinates are the ground state energies of the GDMl 
and GDM2 calculations, respectively, measured from the 
exact one (E'gdmi = “ ^sheii and 

SgdM2 = ggg 

all the three calculations give good ground state ener¬ 
gies: the errors are small relative to the pairing cor¬ 
relation energy. The average values are i?var = 0.045, 
.Egdmi = 0.061, i?GDM 2 = 0.070, and impair = 13.64 
MeV. It is well known that the variational calculation 
finds the best and lowest ground state energy for a set of 
restricted wavefunctions of the form (ED- However, from 
the above average values we see that the energies of the 
GDM wavefunctions and are close to 

the variational minimum {cj)^’^\H\(j>™). 

IV. SUMMARY 

In summary, I derive a physically improved version of 
the GDM pairing theory proposed in Ref. [I|. Both 
versions are checked in a large random ensemble, and 
the accuracy is given statistically in terms of root-mean- 
square derivations from the exact results. Consequently, 
we could consider the theories to be correct and apply 
them with confidence to realistic systems. 

Based on the results, the GDM theories are not as 
accurate as (although close to) the variation principle. 


However, the reduction of computing-time cost is huge 
for large model spaces (by a factor of the number of 
non-degenerate single-particle levels). This should be in¬ 
teresting for ab-initio mean-field calculations of medium 
and heavy nuclei (especially if deformed Nilsson single¬ 
particle levels were used), or the effort in fitting param¬ 
eters of an interaction. 

Comparing the two versions of the GDM theory, we see 
that in general the new one (physically more reasonable) 
is slightly better in calculating pair-transfer amplitudes 
Si, while the old one produces better occupation numbers 
ni and slightly better ground state energy. The reason 
for the latter is still unclear. 

The key approximation of the current GDM methods 
is the “factorization” or “linearization” of the two-body 
density matrix on the pair condensate [Eqs. (|T(I1) and 
(fTTI) ]. It would be interesting to see its validity in other 
circumstances, in particular, whether we could use it 
in the variational formalism when calculating the two- 
body part of the average energy, which would reduce the 
time cost to the same level of the current GDM methods. 
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Michigan State University. 


[1] L. Y. Jia, Phys. Rev. C 88, 044303 (2013). 

[2] A. Bohr, B. R. Mottelson, and D. Pines, Phys. Rev. 110, 
936 (1958). 

[31 S. T. Belyaev, K. Dan. Vidensk. Selsk. Mat. Fys. Medd. 
31, (11) (1959). 

[4] Ricardo A Broglia, and Vladimir Zelevinsky, Fifty Years 
of Nuclear BCS: Pairing in Finite Systems (World Sci¬ 
entific, 2013). 

[5] J. Bardeen, L. N. Cooper, and J. R. Schrieffer, Phys. 
Rev. 106, 162 (1957); Phys. Rev. 108, 1175 (1957). 

[6] A. Bohr and B. Mottelson, Nuclear Structure (Benjamin, 
New York, 1975). 

[7] P. Ring and P. Schuck, The Nuclear Many-Body Problem 
(Springer-Verlag, Berlin, 1980). 

[8] K. Dietrich, H. J. Mang, and J. H. Pradal, Phys. Rev. 
135, B22 (1964). 

[9] L.M. Robledo, G.F. Bertsch, arXiv:1205.4443 [nucl-th] 

( 2012 ). 

[10] H. Flocard and N. Onishi, Annals of Physics 254, 275 
(1997). 

[11] Fabian Braun and Jan von Delft, Phys. Rev. Lett. 81, 
4712 (1998). 

[12] J. Dukelsky and G. Sierra, Phys. Rev. B 61, 12302 

( 2000 ). 

[13] Jan von Delft and D.C. Ralph, Physics Reports 345, 61 

( 2001 ). 

[14] A. Kerman and A. Klein, Phys. Rev. 132, 1326 (1963). 

[15] S.T. Belyaev and V.G. Zelevinsky, Yad. Fiz. 11, 741 


(1970) [Sov. J. Nucl. Phys. 11, 416 (1970)]; Yad. Fiz. 16, 
1195 (1972) [Sov. J. Nucl. Phys. 16, 657 (1973)]; Yad. 
Fiz. 17, 525 (1973) [Sov. J. Nucl. Phys. 17, 269 (1973)]. 

[16] V.G. Zelevinsky, Prog. Theor. Phys. Suppl. 74-75, 251 
(1983). 

[17] M.I. Shtokman, Yad. Fiz. 22, 479 (1975) [Sov. J. Nucl. 
Phys. 22, 247 (1976)]. 

[18] L. Y. Jia, Phys. Rev. C 84, 024318 (2011). 

[19] L. Y. Jia, and V. G. Zelevinsky , Phys. Rev. G 84, 064311 

( 2011 ). 

[20] L. Y. Jia, and V. G. Zelevinsky , Phys. Rev. G 86, 014315 

( 2012 ). 

[21] W.A. Richter, M.G. Van Der Merwe, R.E. Julies, and 
B.A. Brown, Nucl. Phys. A523, 325 (1991). 

[22] H. Molique and J. Dudek, Phys. Rev. G 56, 1795 (1997). 

[23] V. Zelevinsky and A. Volya, Nucl. Phys. A752, 325 
(2005). 

[24] N. Cerf and O. Martin, Phys. Rev. G 47, 2610 (1993); 

[25] N. Cerf, Nucl. Phys. A564, 383 (1993). 

[26] Abhishek Mukherjee, Y. Alhassid, and G. F. Bertsch, 
Phys. Rev. C 83, 014319 (2011). 

[27] R. W. Richardson, Phys. Lett. 3, 277 (1963); Phys. Rev. 
141, 949 (1966). 

[28] J. Dukelsky, S. Pittel, and G. Sierra, Rev. Mod. Phys. 
76, 643 (2004). 

[29] J. Dukelsky, S. Pittel, arXiv:1204.2950 [nucl-th] (2012). 



5 


occupation numbers ni 




FIG. 1: (Color online) Occupation numbers ni of all the 3154 
single-particle levels in the ensemble by different calculations. 
In panel (a), the single-particle levels and the points are in 
one-to-one correspondence, with the horizontal coordinate of 
the point being the exact n\ and the vertical coordinate be¬ 
ing the GDMl m. Similarly for panels (b) and (c), but the 
vertical coordinates are ni of the GDM2 and the variational 
calculations, respectively. See text for details. 
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pair-transfer amplitudes si 



exact exact 


FIG. 2: (Color online) Pair-transfer amplitudes si of all the 3154 single-particle levels in the ensemble by different calculations. 
The points are plotted in the same way as those in Fig. o, but for the transfer amplitudes si. Panels (c) and (d) plot two 
sets of variational calculations (see text for details). 
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single-particle-level groups with different j 
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FIG. 3: (Color online) Results of different calculations 
grouped by the single-particle-level angular momentum j. I 
divide all the 3154 single-particle levels in the ensemble into 
different groups according to their angular momentum j, and 
panel (a) plots the group sizes. Panel (b) plots the root-mean- 
square derivations from the exact results of the occupation 
numbers ni by different calculations (GDMl, GDM2, VAR) 
within each j group. Similarly, panel (c) plots the root-mean- 
square derivations of the pair-transfer amplitudes si by four 
sets of calculations (GDMl, GDM2, VARl, VAR2). See text 
for details. 


















sub-ensembles with different N 
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FIG. 4: (Color online) Results of different calculations in sub- 
ensembles divided by the particle number 2N. The 1000 ex¬ 
amples in the ensemble are divided into sub-ensembles accord¬ 
ing to their particle number 2N, and panel (a) plots the sizes 
of the sub-ensembles. Panels (b) and (c) plot the root-mean- 
square derivations by different calculations of the occupation 
numbers ni and pair-transfer amplitudes si within each N 
sub-ensemble. See text for details. 
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sub-ensembles with different pairing strength 





FIG. 5: (Color online) Results of different calculations in sub- 
ensembles divided by the pairing strength Gmax- The 1000 
examples in the ensemble are divided into 20 sub-ensembles 
according to their pairing strength Gmax, and panel (a) plots 
the sizes of the sub-ensembles. Panels (b) and (c) plot the 
root-mean-square derivations by different calculations of the 
occupation numbers m and pair-transfer amplitudes si within 
each sub-ensemble. See text for details. 
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ground state energies 
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FIG. 6: (Color online) Gronnd state energies of all the 1000 
examples by different calcnlations. In panel (a), the exam¬ 
ples and the points are in one-to-one correspondence, with 
the horizontal coordinate of the point being the pairing cor¬ 
relation energy Bpair, and the vertical coordinate being the 
variational ground state energy iSvar measnred from the ex¬ 
act gronnd state energy. Similarly for panels (b) and (c), 
but the vertical coordinates are ground state energies by the 
GDMl and GDM2 calculations, respectively, measured from 
the exact ground state energy. See text for details. 

















